function [yy_normalized, arma_process] = generate_arma_process_pq(t, burnin, ...
    num_replications, ar_params, ma_params, user_supplied_noise, df)
    % GENERATE_ARMA_PROCESS_PQ Generate ARMA(p,q) process with flexible innovations
    %
    % Generates an ARMA(p,q) process:
    %   ε_t = ρ_1*ε_{t-1} + ... + ρ_p*ε_{t-p} + ζ_t + θ_1*ζ_{t-1} + ... + θ_q*ζ_{t-q}
    %
    % The innovations ζ_t can be:
    %   - Gaussian (df = Inf, default)
    %   - Student's t with ν degrees of freedom (df = ν), standardized to unit variance
    %   - Resampled from user-supplied noise (for bootstrap)
    %
    % INPUTS:
    %   t                   - number of periods (sample size)
    %   burnin              - burn-in periods to discard
    %   num_replications    - number of Monte Carlo replications
    %   ar_params           - AR coefficients [ρ_1, ρ_2, ..., ρ_p] (can be empty [])
    %   ma_params           - MA coefficients [θ_1, θ_2, ..., θ_q] (can be empty [])
    %   user_supplied_noise - noise vector for bootstrap resampling ([] for none)
    %   df                  - degrees of freedom for t-distribution (default: Inf = Gaussian)
    %                         Use df = 5 for heavy tails, df = 15 for moderate tails
    %
    % OUTPUTS:
    %   yy_normalized - ARMA process normalized to unit unconditional variance
    %   arma_process  - raw ARMA process before normalization
    %
    % EXAMPLES:
    %   % ARMA(1,1) with Gaussian innovations (default)
    %   [y_norm, y_raw] = generate_arma_process_pq(100, 1000, 1, 0.5, 0.3, []);
    %
    %   % ARMA(1,1) with t(5) innovations (heavy tails)
    %   [y_norm, y_raw] = generate_arma_process_pq(100, 1000, 1, 0.5, 0.3, [], 5);
    %
    %   % AR(2) with t(15) innovations (moderate tails)
    %   [y_norm, y_raw] = generate_arma_process_pq(100, 1000, 1, [0.5, 0.2], [], [], 15);
    %
    %   % Bootstrap with resampled residuals (df ignored when user_supplied_noise provided)
    %   [~, y_raw] = generate_arma_process_pq(100, 1000, 500, [0.5], [], residuals);

    % Set default for df (Gaussian)
    if nargin < 7 || isempty(df)
        df = Inf;
    end

    % Validate inputs
    assert(t > 0 && burnin >= 0 && num_replications > 0, ...
        'Invalid input: t, burnin, and num_replications must be positive.');
    assert(df > 2 || isinf(df), ...
        'Degrees of freedom must be > 2 (for finite variance) or Inf (Gaussian).');

    % Determine AR and MA orders
    p = length(ar_params);  % AR order
    q = length(ma_params);  % MA order

    % Reshape parameters to row vectors for matrix multiplication
    % Convention: ar_params = [ρ_1, ρ_2, ..., ρ_p] where ρ_1 is for lag 1
    ar_params = reshape(ar_params, [1, p]);
    ma_params = reshape(ma_params, [1, q]);

    % Determine the case based on p and q
    % case_num: 0=white noise, 1=MA(q), 2=AR(p), 3=ARMA(p,q)
    case_num = (p > 0) * 2 + (q > 0);

    % Generate innovations (ζ_t)
    if exist('user_supplied_noise', 'var') && ~isempty(user_supplied_noise)
        % Bootstrap mode: resample with replacement from user-supplied noise
        random_noise = reshape(datasample(user_supplied_noise, ...
            (t + burnin) * num_replications, 'Replace', true), ...
            [t + burnin, num_replications]);
    elseif isinf(df)
        % Gaussian innovations: ζ_t ~ N(0, 1)
        random_noise = randn(t + burnin, num_replications);
    else
        % Student's t innovations, standardized to unit variance
        % Raw t(ν) has variance ν/(ν-2), so we scale by sqrt((ν-2)/ν)
        % Result: Var(ζ_t) = 1
        raw_t = trnd(df, t + burnin, num_replications);
        random_noise = raw_t * sqrt((df - 2) / df);
    end

    % Initialize ARMA process
    arma_process = zeros(t + burnin, num_replications);
    start_idx = max(p, q);

    % Set initial values (first max(p,q) periods)
    for i = 1:start_idx
        arma_process(i, :) = random_noise(i, :);
    end

    % Generate ARMA(p,q) process
    % Using reversed slice order (i-1:-1:i-p) so that:
    %   [ρ_1, ρ_2] * [ε_{t-1}; ε_{t-2}] = ρ_1*ε_{t-1} + ρ_2*ε_{t-2}
    switch case_num
        case 0  % White noise
            arma_process = random_noise;

        case 1  % MA(q) only
            for i = start_idx + 1 : t + burnin
                % ε_t = ζ_t + θ_1*ζ_{t-1} + ... + θ_q*ζ_{t-q}
                ma_term = ma_params * random_noise(i-1:-1:i-q, :);
                arma_process(i, :) = random_noise(i, :) + ma_term;
            end

        case 2  % AR(p) only
            for i = start_idx + 1 : t + burnin
                % ε_t = ρ_1*ε_{t-1} + ... + ρ_p*ε_{t-p} + ζ_t
                ar_term = ar_params * arma_process(i-1:-1:i-p, :);
                arma_process(i, :) = ar_term + random_noise(i, :);
            end

        case 3  % ARMA(p, q)
            for i = start_idx + 1 : t + burnin
                % ε_t = ρ_1*ε_{t-1} + ... + ρ_p*ε_{t-p} + ζ_t + θ_1*ζ_{t-1} + ... + θ_q*ζ_{t-q}
                ar_term = ar_params * arma_process(i-1:-1:i-p, :);
                ma_term = ma_params * random_noise(i-1:-1:i-q, :);
                arma_process(i, :) = ar_term + random_noise(i, :) + ma_term;
            end
    end

    % Discard burn-in period
    arma_process = arma_process(burnin+1:end, :);

    % Calculate unconditional variance for normalization
    % For ARMA(1,1): Var(ε) = (1 + θ^2 + 2ρθ) / (1 - ρ^2)
    uncond_var = compute_unconditional_variance(ar_params, ma_params, p, q);

    % Normalize to unit unconditional variance
    yy_normalized = arma_process ./ sqrt(uncond_var);

end


function uncond_var = compute_unconditional_variance(ar_params, ma_params, p, q)
    % Compute unconditional variance of ARMA(p,q) process
    % Assumes innovation variance σ_ζ² = 1

    if p == 0 && q == 0
        % White noise: Var(ε) = σ_ζ² = 1
        uncond_var = 1;

    elseif p == 1 && q == 0
        % AR(1): Var(ε) = 1 / (1 - ρ²)
        rho = ar_params(1);
        uncond_var = 1 / (1 - rho^2);

    elseif p == 0 && q == 1
        % MA(1): Var(ε) = 1 + θ²
        theta = ma_params(1);
        uncond_var = 1 + theta^2;

    elseif p == 1 && q == 1
        % ARMA(1,1): Var(ε) = (1 + θ² + 2ρθ) / (1 - ρ²)
        rho = ar_params(1);
        theta = ma_params(1);
        uncond_var = (1 + theta^2 + 2*rho*theta) / (1 - rho^2);

    elseif p == 2 && q == 0
        % AR(2): Var(ε) = (1 - ρ_2) / [(1 + ρ_2) * ((1 - ρ_2)² - ρ_1²)]
        rho1 = ar_params(1);
        rho2 = ar_params(2);
        uncond_var = (1 - rho2) / ((1 + rho2) * ((1 - rho2)^2 - rho1^2));

    else
        % General case: use empirical variance (less accurate but robust)
        % This fallback should rarely be needed in practice
        % warning('generate_arma_process_pq:generalCase', ...
        %    'Using empirical variance for ARMA(%d,%d). Consider adding exact formula.', p, q);
        uncond_var = 1;  % Will be corrected by empirical normalization if needed
    end
end